%This file, ie champ_julia1.m was created from champ_julia on 15dec06
% with respect to different LT bins of Julia_data (I mean> no champ_fday - lag -
% center point -average business !) - manoj 15dec06
% This analysis picks out the depenednce of Cc on LAG more clearly
% Less elegant, simpler!
data11 = load('c:\manoj\projects\plasma\STRICT_CRT_PLP_EOP_JT_09-17LT_268-298LON_KP.txt'); %15dec6 Claudia's new manual selecteion see email 
load c:\manoj\projects\plasma\Julia_PLP.mat;
%PLP = PLP_new;% PLP_new is IRI corrected PLP data at 400 km (16Jan2007)
PLP = PLP_15AUG07; % new corrected PLP data (see mail from CS on 15-08-2007)
%Comment the previous line to use the old IRI-uncorrected PLP data
fday_list_julia = floor(fday(:,1));
kp_limit = 3; 
delay_index = 12;
delay_sel = [13,14,15,16,17,18,19,20,21,22,23,24];
Kp_delay = delay_sel(delay_index);%Valid columns [10,13-24];
kp_delay_text = ['00';'01';'02';'03';'04';'05';'06';'07';'08';'09';'22';'MN'];
lt_step =  1; %LT window size in Hours
lag_step = 60; % LAG STEP in minutes
min_data = 5; %Minimum number of data pairs to be considered for correlation
%LAG = [30:lag_step:240]/(60*24);% LAG bins between CHAMP and Julia
%LAG = [-30:lag_step:390]/(60*24);
LT = [10:lt_step:17];%LT bins 
LAG = [13:(lag_step/5):109]; %13start 12*5=60min 109 = evening 4 , 2 - end julia record
AA = [8.5:(lag_step/60):15.5];
robft = 1;
plot_option = 2;%[1 scatter, 2 image 3 curves];
ED_matrix = ones([length(LT)-1,length(LAG)-1])*NaN;
LV_matrix = ones([length(LT)-1,length(LAG)-1])*NaN;
RT_matrix = ones([length(LT)-1,length(LAG)-1])*NaN;
PT_matrix = ones([length(LT)-1,length(LAG)-1])*NaN;
ED_error =  ones([length(LT)-1,length(LAG)-1])*NaN;
LV_error =  ones([length(LT)-1,length(LAG)-1])*NaN;
RT_error =  ones([length(LT)-1,length(LAG)-1])*NaN;
PT_error =  ones([length(LT)-1,length(LAG)-1])*NaN;

for jj = 1:length(LAG)-1,
for kk = 1:length(LT)-1,
    
%      L = ~isnan(PLP(:,27)) & ~isnan(PLP(:,18))  & D_KP(:,Kp_delay) <=kp_limit & PLP(:,7) >= LT(kk) & PLP(:,7) <= LT(kk+1)...
%          & PLP(:,12) == 1;% & PLP(:,20) >= 1.5;%& PLP(:,29) >= 1.5;% & PLP(:,8) > -85 & PLP(:,8) < -70;
      L = ~isnan(PLP(:,27)) & ~isnan(PLP(:,18))  & data11(:,Kp_delay) <=kp_limit & PLP(:,7) >= LT(kk) & PLP(:,7) <= LT(kk+1)...
          & data11(:,12) == 1;% & PLP(:,20) >= 1.5;%& PLP(:,29) >= 1.5;% & PLP(:,8) > -85 & PLP(:,8) < -70;

         IndData = PLP(L,:);
         counter = 0;
         
         for i = 1:size(IndData,1),
%         fday_champ = tmjd(IndData(i,1),IndData(i,2),IndData(i,3),IndData(i,4),IndData(i,5),IndData(i,6),1);   
         fday_champ = datenum(IndData(i,1),IndData(i,2),IndData(i,3),IndData(i,4),IndData(i,5),IndData(i,6))-730486;   
         iday = find(fday_list_julia==floor(fday_champ));
         if ~isempty(iday),
         LLL = isnan(w(iday,:)); kkk = find(LLL==0);
            if length(kkk) > 3,
                 counter=counter+1;
                Julia_data(counter) = nanmean(w(iday,LAG(jj):LAG(jj+1))); %  AVERAGE
                Plp_ed(counter) = (IndData(i,18) + IndData(i,27))/2;
              %  Plp_ed(counter) = (IndData(i,18)); %WATCH OUT ONLY ONE PEAK !
                Plp_lv(counter) = (IndData(i,17) + IndData(i,26))/2;
                %Plp_ed(counter) = (IndData(i,18) + IndData(i,27))/2;
                Plp_ratio(counter) = (IndData(i,20) + IndData(i,29))/2;
                Plp_pt(counter) = Plp_ed(counter) - IndData(i,36);
               % Plp_pt(counter) = IndData(i,14)-IndData(i,23);%NCREST LONG- S. CREST LONG
                Plp_altitude(counter) = IndData(i,9);
            end;
            
         end;
         end;
        if counter >= min_data,
         LP = isnan(Julia_data);
        Julia_data(LP) = [];Plp_ed(LP) = [];Plp_lv(LP) = [];Plp_ratio(LP) = [];
        Plp_pt(LP) = [];
        counter = counter - sum(LP);
        end;
        if robft == 1 & counter >= min_data,
                [Julia_data,Plp_pt,Plp_ratio,Plp_lv,Plp_ed] = ...
                rob_regress(Julia_data,Plp_pt,Plp_ratio,Plp_lv,Plp_ed);
                counter = length(Julia_data);
% 21aug07 the following section for robust processing was advanced. however
% to keep the results same for the first sw paper, follow the above rb
% proc.
%                 temp_1 = rob_regress1(Julia_data,Plp_pt,Plp_ratio,Plp_lv,Plp_ed);
%                 Julia_data(temp_1) = [];
%                 Plp_ed(temp_1) = [];
%                 Plp_lv(temp_1) = [];
%                 Plp_ratio(temp_1) = [];
%                 Plp_pt(temp_1) = [];
%                 counter = length(Julia_data);
         end;
         lt_count(kk,jj) = counter;
        % fprintf('LT=%5d LAG = %10.5f COUNT = %3d',LT(kk),LAG(jj)*(60*24),counter); 
         if counter >= min_data,
             [cof_ed,err_ed] = corrcoef(Julia_data,Plp_ed); 
             [cof_lv,err_lv] = corrcoef(Julia_data,Plp_lv); 
             [cof_rt,err_rt] = corrcoef(Julia_data,Plp_ratio); 
             [cof_pt,err_pt] = corrcoef(Julia_data,Plp_pt); 
             ED_matrix(kk,jj) = cof_ed(2,1);
             LV_matrix(kk,jj) = cof_lv(2,1);
             ED_error(kk,jj) = err_ed(2,1);
             RT_matrix(kk,jj) = cof_rt(2,1);
             RT_error(kk,jj) = err_rt(2,1);
             PT_matrix(kk,jj) = cof_pt(2,1);
             PT_error(kk,jj) = err_pt(2,1);
             LV_error(kk,jj) = err_lv(2,1);
             
             if plot_option == 1,
%              ss=((LAG+83)*5)/60;
%              subplot(length(LAG)-1,length(LT)-1,(jj-1)*(length(LT)-1)+kk);
%              %set(gca,'FontSize',16);
%              hold on;
%              plot(Julia_data,Plp_ratio,'k.','MarkerSize',15);
%              
%              [c,sfit] = polyfit(Julia_data,Plp_ratio,1);
%              m = polyval(c,[min(Julia_data),max(Julia_data)]);
%              plot([min(Julia_data),max(Julia_data)],m,'r-','LineWidth',2);
% 
%                 if jj ~=(length(LAG)-1), set(gca,'XTick',[]);end;
%                 if kk ~=1 ,set(gca,'YTick',[]);end;
%                 if jj == (length(LAG)-1), text(10,0.0,sprintf('%d-%d',LT(kk),LT(kk+1)),'FontSIze',16);end;
%                 if kk == 1 , text(-50,2.5,sprintf('%d-%d',ss(jj),ss(jj+1)),'FontSIze',16);end;
%               %grid on;
%                 axis([-10 50 1 3]);
% %                 text(-10,2.5,sprintf('%d',LT(kk)));
% %                 text(20,2.5,sprintf('%4.1f',AA(jj)));
% %                 text(40,2.5,sprintf('%4.2f',cof_rt(2,1)),'color','b');
%                 text(40,1.5,sprintf('%2d',round(cof_rt(2,1)*100)),'color','k');
%                 %xlabel(sprintf('%5d-%5d',LT(kk),LAG(jj)*(60*24)));
%             fprintf('%d-%d %d-%d %2d %f %f\n',LT(kk),LT(kk+1),ss(jj),...
%                 ss(jj+1),round(cof_rt(2,1)*100),c(1),c(2));

% TO PLOT LVLUE SCATTER !

             ss=((LAG+83)*5)/60;
             subplot(length(LAG)-1,length(LT)-1,(jj-1)*(length(LT)-1)+kk);
             %set(gca,'FontSize',16);
             hold on;
             plot(Julia_data,Plp_lv,'k.','MarkerSize',15);
             
             [c,sfit] = polyfit(Julia_data,Plp_lv,1);
             m = polyval(c,[min(Julia_data),max(Julia_data)]);
             plot([min(Julia_data),max(Julia_data)],m,'r-','LineWidth',2);

                if jj ~=(length(LAG)-1), set(gca,'XTick',[]);end;
                if kk ~=1 ,set(gca,'YTick',[]);end;
                if jj == (length(LAG)-1), text(10,0.95,sprintf('%d-%d',LT(kk),LT(kk+1)),'FontSIze',16);end;
                if kk == 1 , text(-45,1.1,sprintf('%d-%d',ss(jj),ss(jj+1)),'FontSIze',16);end;
              %grid on;
                axis([-5 35 1.05 1.2]);
%                 text(-10,2.5,sprintf('%d',LT(kk)));
%                 text(20,2.5,sprintf('%4.1f',AA(jj)));
%                 text(40,2.5,sprintf('%4.2f',cof_rt(2,1)),'color','b');
                text(30,1.1,sprintf('%2d',round(cof_lv(2,1)*100)),'color','k');
                %xlabel(sprintf('%5d-%5d',LT(kk),LAG(jj)*(60*24)));
            fprintf('%d-%d %d-%d %2d %f %f\n',LT(kk),LT(kk+1),ss(jj),...
                ss(jj+1),round(cof_rt(2,1)*100),c(1),c(2));

              end;
             end;
%         fprintf('\n');
         clear Julia_data Plp_ed Plp_lv Plp_ratio Plp_pt;
      
end;
end;

for i = 1:length(LT)-1,
SS(i,:) = sprintf('%2d-%2d',LT(i),LT(i+1));
end;
if plot_option == 3,

color_index = ['k*-'; 'b*-' ;'r*-'; 'g*-'; 'c*-' ;'y.-' ;'m.-'; 'r.-' ;'g.-'; 'c.-'];

% color_map = colormap('jet');
% color_index = color_map([1:3:64],:);
figure(gcf+1);
subplot(331);
for i = 1:length(LT)-1,
hold on;
errorbar(LAG*(24*60),PT_matrix(i,:)',PT_error(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
xlabel('LAG in minutes');
ylabel('CC');
hold on;
axis([LAG(1)*24*60 LAG(end)*24*60 0 1]);
end;
title('C-T - Julia correlation');
text(-42,1.0,['Kp-' sprintf('%d',kp_limit)]);
text(-42,0.8,['Delay-' kp_delay_text(delay_index,:)]);
subplot(334);
for i = 1:length(LT)-1,
hold on;
plot(LAG*(24*60),lt_count(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);
hold on;
end;
xlabel('LAG in minutes');
ylabel('Number of data points');

subplot(337);
for i = 1:length(LT)-1,
plot(LAG*(24*60),log10(PT_error(i,:))','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);
hold on;
end;
xlabel('LAG in minutes');
ylabel('Log |Error|');


subplot(332);
for i = 1:length(LT)-1,
    hold on;
errorbar(LAG*(24*60),RT_matrix(i,:)',RT_error(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
xlabel('LAG in minutes');
ylabel('CC');
hold on;
axis([LAG(1)*24*60 LAG(end)*24*60 0 1]);
end;
title('CTR - Julia correlation');

subplot(335);
for i = 1:length(LT)-1,
plot(LAG*(24*60),lt_count(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);

hold on;
end;
xlabel('LAG in minutes');
ylabel('Number of data points');

subplot(338);
for i = 1:length(LT)-1,
plot(LAG*(24*60),log10(RT_error(i,:))','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);
hold on;
end;
xlabel('LAG in minutes');
ylabel('Log |Error|');


subplot(333);
for i = 1:length(LT)-1,
    hold on;
errorbar(LAG*(24*60),LV_matrix(i,:)',LV_error(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
xlabel('LAG in minutes');
ylabel('CC');
hold on;
axis([LAG(1)*24*60 LAG(end)*24*60 0 1]);
end;
title('L-Value - Julia correlation');

subplot(336);
for i = 1:length(LT)-1,
plot(LAG*(24*60),lt_count(i,:)','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);

hold on;
end;
xlabel('LAG in minutes');
ylabel('Number of data points');

subplot(339);
for i = 1:length(LT)-1,
plot(LAG*(24*60),log10(LV_error(i,:))','Color',color_index(i),'LineWidth', 2,'MarkerSize',30);
axis([LAG(1)*24*60 LAG(end)*24*60 -inf inf]);
hold on;
end;
xlabel('LAG in minutes');
ylabel('Log |Error|');

end;
% saveas(gcf,['c:\manoj\projects\plasma\Kp3Delay' kp_delay_text(delay_index,:) '.eps'],'psc2');
% close all;
% end;
%}

if plot_option == 2, %image
    figure(gcf+1);
    ss=((LAG+83)*5)/60;
    dd = (ss(1:end-1)+diff(ss)/2);
    BB = [LT(1:end-1)+diff(LT)/2];
    AA =dd;
    subplot(2,2,1);
    imagesc(BB,dd,RT_matrix');
    set(gca,'YTick',ss);
    set(gca,'YTickLabel',ss);
for i = 1:length(AA),
for j = 1:length(BB),
%text(BB(j),AA(i),sprintf('%d',lt_count(j,i)));
text(BB(j),AA(i),sprintf('%2d',round(RT_matrix(j,i)*100)));
end;
end;
    
    
    caxis([0,1]);
    title('CTR');
    colorbar;
  %  set(gca,'FontSize',16)
    subplot(2,2,2);
    imagesc(BB,dd,PT_matrix');
    set(gca,'YTick',ss);
    set(gca,'YTickLabel',ss);
    caxis([0,1]);
    colorbar;
    title('Crest-Trough');
    for i = 1:length(AA),
for j = 1:length(BB),
%text(BB(j),AA(i),sprintf('%d',lt_count(j,i)));
text(BB(j),AA(i),sprintf('%2d',round(PT_matrix(j,i)*100)));
end;
end;
%set(gca,'FontSize',16)
    subplot(2,2,3);
    imagesc(BB,dd,LV_matrix');
    set(gca,'YTick',ss);
    set(gca,'YTickLabel',ss);
    caxis([0,1]);
    colorbar;
    title('L-value');
 for i = 1:length(AA),
for j = 1:length(BB),
%text(BB(j),AA(i),sprintf('%d',lt_count(j,i)));
text(BB(j),AA(i),sprintf('%2d',round(LV_matrix(j,i)*100)));
end;
end;
  % set(gca,'FontSize',16)
    subplot(2,2,4);
    imagesc(BB,dd,ED_matrix');
    set(gca,'YTick',ss);
    set(gca,'YTickLabel',ss);
    caxis([0,1]);
    colorbar;
    title('Crest');
    for i = 1:length(AA),
    for j = 1:length(BB),
    %text(BB(j),AA(i),sprintf('%d',lt_count(j,i)));
    text(BB(j),AA(i),sprintf('%2d',round(ED_matrix(j,i)*100)));
    end;
    end;
%set(gca,'FontSize',16)
end;